#ifndef _TemplateElement_
#define _TemplateElement_

#include "Dofs.h"
#include "GaussianPoint.h"
#include "Grid.h"
#include "Mesh.h"
#include <math.h>
/**
 * @brief  建立一个参考单元，只考虑二维情况
 */
template <unsigned int DIM>
class TemplateElement
{
public:
    /**
     * @brief 默认构造函数
     *
     */
    TemplateElement() = default;
    /**
     * @brief 默认析构函数
     *
     */
    virtual ~TemplateElement() = default;

    /**
     * @brief 从当前绑定的Grid中获取自由度信息
     *
     * @param i 指定序号；
     * @return Dofs<DIM>&　自由度的引用
     */
    virtual Dofs<DIM> &getDofsList(int i) = 0;
    /**
     * @brief 从当前绑定的Grid中获取第ｉ个自由度对应的全局编号；
     *
     * @param i 指定序号；
     * @return int 全局编号；
     */
    virtual int getGlobalIndex(int i) = 0;
    /**
     * @brief 返回本单元对应第i个高斯点坐标信息；
     *
     * @param i 指定序号
     * @return Point<DIM>& 高斯点
     */
    virtual Point<DIM> &GaussionPoint(int i) = 0;
    /**
     * @brief 返回本单元对应第i个高斯点权重信息；
     *
     * @param i 指定序号
     * @return double 高斯点对应权重的数值；
     */
    virtual double GaussionWeight(int i) = 0;
    /**
     * @brief 该类型单元对应的基函数；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i　基函数序数
     * @return double 该基函数的值；
     */
    virtual double phi(double xi, double eta, int i) = 0;
    /**
     * @brief 该类型单元对应的基函数对xi的偏导数；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 该偏导数的值
     */
    virtual double phi_xi(double xi, double eta, int i) = 0;
    /**
     * @brief 该类型单元对应的基函数对eta的偏导数；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 该偏导数的值
     */
    virtual double phi_eta(double xi, double eta, int i) = 0;
    /**
     * @brief 任意单元到该模板单元仿射变换对应的雅克比矩阵行列式在(xi,eta)处的值；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @return double 值
     */
    virtual double det_Jacobi(double xi, double eta) = 0;
    /**
     * @brief 基函数对x的偏导数；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 值
     */
    virtual double phi_x(double xi, double eta, int i) = 0;
    /**
     * @brief 基函数对y的偏导数；
     *
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 值
     */
    virtual double phi_y(double xi, double eta, int i) = 0;
    /**
     * @brief 计算并返回第ｉ个基函数在xi，eta点处的梯度向量；
     *
     * @param xi
     * @param eta
     * @param i
     * @return std::vector<double> 梯度向量；
     */
    virtual std::vector<double> gradient(double xi, double eta, int i) = 0;
    /**
     * @brief 标准单元中的点在真实单元的x坐标
     *
     * @param xi 坐标１
     * @param eta 坐标2
     * @return double 值
     */
    virtual double Global_x(double xi, double eta) = 0;
    /**
     * @brief 标准单元中的点在真实单元的y坐标
     *
     * @param xi 坐标１
     * @param eta 坐标2
     * @return double 值
     */
    virtual double Global_y(double xi, double eta) = 0;
    /**
     * @brief 该单元自由度个数；
     *
     * @return int
     */
    virtual int n_Dofs() = 0;
    /**
     * @brief 高斯点个数；
     *
     * @return int
     */
    virtual int n_GaussPnt() = 0;
    /**
     * @brief 计算并返回该单元的面积或体积；
     *
     * @return double
     */
    virtual double volume() = 0;
    /**
     * @brief 获取grid的自由度等信息
     * @details 将grid中的dofs赋给_dofslist
     */
    virtual void reinit(Grid<DIM> &_grid, Mesh<DIM> &_mesh) = 0;
    // virtual void setMesh(Mesh<DIM>& _mesh) = 0;
    virtual vector<Dofs<2>> &getBoundaryDofsonGrid() = 0;

protected:
    virtual double xi_x(double xi, double eta) = 0;
    virtual double xi_y(double xi, double eta) = 0;
    virtual double eta_x(double xi, double eta) = 0;
    virtual double eta_y(double xi, double eta) = 0;
};

/**
 *@brief 一阶四边形模板单元
 */
class Quadrilateral_1_Element : public TemplateElement<2>
{
    typedef std::vector<Dofs<2>> DofsList;
    typedef int DofsNum;

public:
    Quadrilateral_1_Element();
    ~Quadrilateral_1_Element() = default;
    virtual Dofs<2> &getDofsList(int i);
    virtual int getGlobalIndex(int i);
    virtual Point<2> &GaussionPoint(int i) { return GaussInfo.GaussianPointInfo(i); };
    virtual double GaussionWeight(int i) { return GaussInfo.GaussianWeightInfo(i); };
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i);
    virtual double phi_eta(double xi, double eta, int i);
    virtual double det_Jacobi(double xi, double eta);
    virtual double xi_x(double xi, double eta);
    virtual double xi_y(double xi, double eta);
    virtual double eta_x(double xi, double eta);
    virtual double eta_y(double xi, double eta);
    virtual double phi_x(double xi, double eta, int i);
    virtual double phi_y(double xi, double eta, int i);
    virtual std::vector<double> gradient(double xi, double eta, int i);
    virtual double Global_x(double xi, double eta);
    virtual double Global_y(double xi, double eta);
    virtual double volume() { return 4 * det_Jacobi(0, 0); };
    virtual int n_Dofs() { return 4; };
    virtual int n_GaussPnt() { return GaussInfo.GaussianPointNumInfo(); };
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh);
    // virtual void setMesh(Mesh<2>& _mesh);
    virtual vector<Dofs<2>> &getBoundaryDofsonGrid();

private:
    DofsList dofslist;                            /*< 存储绑定的Grid中自由度信息*/
    DofsList boundarydofslist;                    /*< 存储绑定的Grid中边界自由度信息*/
    Quadrilateral_Element_GaussianInfo GaussInfo; /*<Gauss积分所需：积分点和积分权重*/
};

// void Quadrilateral_1_Element::setMesh(Mesh<2>& _mesh)
// {
//     mesh.getSegment() = _mesh.getSegment();
//     mesh.getGrid() = _mesh.getGrid();
// }

Quadrilateral_1_Element::Quadrilateral_1_Element()
{
    dofslist.resize(4);
    GaussInfo.LoadGaussianInfo(3); // two dimensions,tensor-product rules
}

Dofs<2> &Quadrilateral_1_Element::getDofsList(int i)
{
    return dofslist[i];
}

double Quadrilateral_1_Element::phi(double xi, double eta, int i)
{
    switch (i)
    {
    case 0:
        return (xi - 1) * (eta - 1) / 4.0;
    case 1:
        return -(xi + 1) * (eta - 1) / 4.0;
    case 2:
        return (xi + 1) * (eta + 1) / 4.0;
    case 3:
        return -(xi - 1) * (eta + 1) / 4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::phi_xi(double xi, double eta, int i)
{
    switch (i)
    {
    case 0:
        return (eta - 1) / 4.0;
    case 1:
        return -(eta - 1) / 4.0;
    case 2:
        return (eta + 1) / 4.0;
    case 3:
        return -(eta + 1) / 4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::phi_eta(double xi, double eta, int i)
{
    switch (i)
    {
    case 0:
        return (xi - 1) / 4.0;
    case 1:
        return -(xi + 1) / 4.0;
    case 2:
        return (xi + 1) / 4.0;
    case 3:
        return -(xi - 1) / 4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::det_Jacobi(double xi, double eta)
{
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 0; i < 4; i++)
    {
        J11 = J11 + dofslist[i][0] * phi_xi(xi, eta, i);
        J12 = J12 + dofslist[i][1] * phi_xi(xi, eta, i);
        J21 = J21 + dofslist[i][0] * phi_eta(xi, eta, i);
        J22 = J22 + dofslist[i][1] * phi_eta(xi, eta, i);
    }
    return (J11 * J22 - J12 * J21);
}

double Quadrilateral_1_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 0; i < 4; i++)
        a = a + dofslist[i][1] * phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 0; i < 4; i++)
        a = a + dofslist[i][0] * phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 0; i < 4; i++)
        a = a + dofslist[i][1] * phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::eta_y(double xi, double eta)
{
    double a = 0;
    for (int i = 0; i < 4; i++)
        a = a + dofslist[i][0] * phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::phi_x(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i) * xi_x(xi, eta) + phi_eta(xi, eta, i) * eta_x(xi, eta);
}

double Quadrilateral_1_Element::phi_y(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i) * xi_y(xi, eta) + phi_eta(xi, eta, i) * eta_y(xi, eta);
}

std::vector<double> Quadrilateral_1_Element::gradient(double xi, double eta, int i)
{
    std::vector<double> gradient({phi_x(xi, eta, i), phi_y(xi, eta, i)});
    return gradient;
}
double Quadrilateral_1_Element::Global_x(double xi, double eta)
{
    double a = 0;
    for (int i = 0; i < 4; i++)
    {
        a += dofslist[i][0] * phi(xi, eta, i);
    }
    return a;
}

double Quadrilateral_1_Element::Global_y(double xi, double eta)
{
    double a = 0;
    for (int i = 0; i < 4; i++)
    {
        a += dofslist[i][1] * phi(xi, eta, i);
    }
    return a;
}

int Quadrilateral_1_Element::getGlobalIndex(int i)
{
    return dofslist[i].getGlobalIndex();
}

void Quadrilateral_1_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh)
{
    dofslist.clear();
    vector<Point<2>> p;
    p.resize(4);
    for (int i = 0; i < 4; i++)
        p[i] = _grid.getVertices(i);
    int j = (int)round((p[0][0] + 1) * 0.5 * _mesh.getSegment()[0]);
    int k = (int)round((p[0][1] + 1) * 0.5 * _mesh.getSegment()[1]);
    int i = j + k * _mesh.getSegment()[0];
    // Point<2> p1 = {Global_x(-1, -1), Global_y(-1, -1)};
    // Point<2> p2 = {Global_x(1, -1), Global_y(1, -1)};
    // Point<2> p3 = {Global_x(1, 1), Global_y(1, 1)};
    // Point<2> p4 = {Global_x(-1, 1), Global_y(-1, 1)};
    for (int k = 0; k < 4; k++)
    {
        Dofs<2> q = {p[k], _mesh.getIndexofGrid(i)[k].first};
        dofslist.push_back(q);
        if (_mesh.getIndexofGrid(i)[k].second != 0)
            boundarydofslist.push_back(q);
    }
};

vector<Dofs<2>> &Quadrilateral_1_Element::getBoundaryDofsonGrid()
{
    return boundarydofslist;
}

class Quadrilateral_2_Element : public TemplateElement<2>
{
    typedef std::vector<Dofs<2>> DofsList;
    typedef int DofsNum;

public:
    Quadrilateral_2_Element();
    ~Quadrilateral_2_Element(){};
    virtual Dofs<2> &getDofsList(int i);
    virtual int getGlobalIndex(int i);
    virtual Point<2> &GaussionPoint(int i) { return GaussInfo.GaussianPointInfo(i); };
    virtual double GaussionWeight(int i) { return GaussInfo.GaussianWeightInfo(i); };
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i);
    virtual double phi_eta(double xi, double eta, int i);
    virtual double det_Jacobi(double xi, double eta);
    virtual double xi_x(double xi, double eta);
    virtual double xi_y(double xi, double eta);
    virtual double eta_x(double xi, double eta);
    virtual double eta_y(double xi, double eta);
    virtual double phi_x(double xi, double eta, int i);
    virtual double phi_y(double xi, double eta, int i);
    virtual std::vector<double> gradient(double xi, double eta, int i);
    virtual double volume() { return 4 * det_Jacobi(0, 0); };
    virtual double Global_x(double xi, double eta);
    virtual double Global_y(double xi, double eta);
    virtual int n_Dofs() { return 9; };
    virtual int n_GaussPnt() { return GaussInfo.GaussianPointNumInfo(); };
    virtual void reinit(Grid<2> &_grid, Mesh<2> &_mesh);
    virtual vector<Dofs<2>> &getBoundaryDofsonGrid();

private:
    DofsList dofslist;                            /*< 存储绑定的Grid中自由度信息*/
    DofsList boundarydofslist;                    /*< 存储绑定的Grid中边界自由度信息*/
    Quadrilateral_Element_GaussianInfo GaussInfo; /*<Gauss积分所需：积分点和积分权重*/
};

Quadrilateral_2_Element::Quadrilateral_2_Element()
{
    dofslist.resize(9);
    GaussInfo.LoadGaussianInfo(5);
}

Dofs<2> &Quadrilateral_2_Element::getDofsList(int i)
{
    return dofslist[i];
}

double Quadrilateral_2_Element::phi(double xi, double eta, int i)
{
    /// 注意基函数顺序与打点分布的一致性
    switch (i)
    {
    case 0:
        return (1 / 4.0) * xi * eta * (1.0 - xi) * (1.0 - eta);
    case 1:
        return -(1 / 2.0) * (1.0 - eta) * eta * (1.0 - xi * xi);
    case 2:
        return -(1 / 4.0) * xi * eta * (1.0 + xi) * (1.0 - eta);
    case 3:
        return (1 / 2.0) * (1.0 + xi) * xi * (1.0 - eta * eta);
    case 4:
        return (1 / 4.0) * xi * eta * (1.0 + xi) * (1.0 + eta);
    case 5:
        return (1 / 2.0) * (1.0 + eta) * eta * (1.0 - xi * xi);
    case 6:
        return -(1 / 4.0) * xi * eta * (1.0 - xi) * (1.0 + eta);
    case 7:
        return -(1 / 2.0) * (1.0 - xi) * xi * (1.0 - eta * eta);
    case 8:
        return (1.0 - xi * xi) * (1.0 - eta * eta);
    }
    return 0;
}

double Quadrilateral_2_Element::phi_xi(double xi, double eta, int i)
{
    switch (i)
    {
    case 0:
        return (1 / 4.0) * (1.0 - 2.0 * xi) * eta * (1.0 - eta);
    case 1:
        return xi * (1.0 - eta) * eta;
    case 2:
        return -(1 / 4.0) * (1.0 + 2.0 * xi) * eta * (1.0 - eta);
    case 3:
        return (1 / 2.0) * (1.0 + 2.0 * xi) * (1.0 - eta * eta);
    case 4:
        return (1 / 4.0) * (1.0 + 2.0 * xi) * eta * (1.0 + eta);
    case 5:
        return -xi * (1.0 + eta) * eta;
    case 6:
        return -(1 / 4.0) * (1.0 - 2.0 * xi) * eta * (1.0 + eta);
    case 7:
        return -(1 / 2.0) * (1.0 - 2.0 * xi) * (1.0 - eta * eta);
    case 8:
        return -2.0 * xi * (1.0 - eta * eta);
    }
    return 0;
}

double Quadrilateral_2_Element::phi_eta(double xi, double eta, int i)
{
    switch (i)
    {
    case 0:
        return (1 / 4.0) * xi * (1.0 - xi) * (1.0 - 2.0 * eta);
    case 1:
        return -(1 / 2.0) * (1.0 - xi * xi) * (1.0 - 2.0 * eta);
    case 2:
        return -(1 / 4.0) * xi * (1.0 + xi) * (1 - 2.0 * eta);
    case 3:
        return -xi * (1.0 + xi) * eta;
    case 4:
        return (1 / 4.0) * xi * (1.0 + xi) * (1.0 + 2.0 * eta);
    case 5:
        return (1 / 2.0) * (1.0 - xi * xi) * (1.0 + 2.0 * eta);
    case 6:
        return -(1 / 4.0) * xi * (1.0 - xi) * (1 + 2.0 * eta);
    case 7:
        return xi * (1.0 - xi) * eta;
    case 8:
        return -2.0 * (1.0 - xi * xi) * eta;
    }
    return 0;
}

double Quadrilateral_2_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 0; i < 9; i++)
        a = a + dofslist[i][1] * phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 0; i < 9; i++)
        a = a + dofslist[i][0] * phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}
double Quadrilateral_2_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 0; i < 9; i++)
        a = a + dofslist[i][1] * phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::eta_y(double xi, double eta)
{
    double a = 0;
    for (int i = 0; i < 9; i++)
        a = a + dofslist[i][0] * phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::phi_x(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i) * xi_x(xi, eta) + phi_eta(xi, eta, i) * eta_x(xi, eta);
}

double Quadrilateral_2_Element::phi_y(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i) * xi_y(xi, eta) + phi_eta(xi, eta, i) * eta_y(xi, eta);
}

std::vector<double> Quadrilateral_2_Element::gradient(double xi, double eta, int i)
{
    std::vector<double> gradient({phi_x(xi, eta, i), phi_y(xi, eta, i)});
    return gradient;
}

double Quadrilateral_2_Element::det_Jacobi(double xi, double eta)
{
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 0; i < 9; i++)
    {
        J11 = J11 + dofslist[i][0] * phi_xi(xi, eta, i);
        J12 = J12 + dofslist[i][1] * phi_xi(xi, eta, i);
        J21 = J21 + dofslist[i][0] * phi_eta(xi, eta, i);
        J22 = J22 + dofslist[i][1] * phi_eta(xi, eta, i);
    }
    return (J11 * J22 - J12 * J21);
}

double Quadrilateral_2_Element::Global_x(double xi, double eta)
{
    double a = 0;
    a += dofslist[0][0] * 0.25 * (1 - xi) * (1 - eta);
    a += dofslist[2][0] * 0.25 * (1 + xi) * (1 - eta);
    a += dofslist[4][0] * 0.25 * (1 + xi) * (1 + eta);
    a += dofslist[6][0] * 0.25 * (1 - xi) * (1 + eta);
    return a;
}

double Quadrilateral_2_Element::Global_y(double xi, double eta)
{
    double a = 0;
    a += dofslist[0][1] * 0.25 * (1 - xi) * (1 - eta);
    a += dofslist[2][1] * 0.25 * (1 + xi) * (1 - eta);
    a += dofslist[4][1] * 0.25 * (1 + xi) * (1 + eta);
    a += dofslist[6][1] * 0.25 * (1 - xi) * (1 + eta);
    return a;
}

int Quadrilateral_2_Element::getGlobalIndex(int i)
{
    return dofslist[i].getGlobalIndex();
}

void Quadrilateral_2_Element::reinit(Grid<2> &_grid, Mesh<2> &_mesh)
{
    dofslist.clear();
    vector<Point<2>> p;
    p.resize(4);
    for (int i = 0; i < 4; i++)
        p[i] = _grid.getVertices(i);
    int j = (int)round((p[0][0] + 1) * 0.5 * _mesh.getSegment()[0]);
    int k = (int)round((p[0][1] + 1) * 0.5 * _mesh.getSegment()[1]);
    int i = j + k * _mesh.getSegment()[0];
    vector<Dofs<2>> q;
    q.resize(9);
    q[0] = {p[0], _mesh.getIndexofGrid(i)[0].first};
    q[1] = {(p[0] + p[1]) / 2, _mesh.getIndexofGrid(i)[1].first};
    q[2] = {p[1], _mesh.getIndexofGrid(i)[2].first};
    q[3] = {(p[1] + p[2]) / 2, _mesh.getIndexofGrid(i)[3].first};
    q[4] = {p[2], _mesh.getIndexofGrid(i)[4].first};
    q[5] = {(p[2] + p[3]) / 2, _mesh.getIndexofGrid(i)[5].first};
    q[6] = {p[3], _mesh.getIndexofGrid(i)[6].first};
    q[7] = {(p[3] + p[0]) / 2, _mesh.getIndexofGrid(i)[7].first};
    q[8] = {(p[0] + p[1] + p[2] + p[3]) / 4, _mesh.getIndexofGrid(i)[8].first};
    for (int k = 0; k < 9; k++)
    {
        dofslist.push_back(q[k]);
        if (_mesh.getIndexofGrid(i)[k].second != 0)
            boundarydofslist.push_back(q[k]);
    }
}

vector<Dofs<2>> &Quadrilateral_2_Element::getBoundaryDofsonGrid()
{
    return boundarydofslist;
}

#else
// Do Nothing.
#endif
